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Abstract 

We tackle the modeling of threshold exceedances in asymptotically independent stochastic 
processes by constructions based on Laplace random fields. Defined as mixtures of Gaussian ran¬ 
dom fields with an exponential variable embedded for the variance, these processes possess useful 
asymptotic properties while remaining statistically convenient. Univariate Laplace distribution 
tails are part of the limiting generalized Pareto distributions for threshold exceedances. After 
normalizing marginal distributions in data, a standard Laplace field can be used to capture spatial 
dependence among extremes. Asymptotic properties of Laplace fields are explored and compared 
to the classical framework of asymptotic dependence. Multivariate joint tail decay rates are slower 
than for Gaussian fields with the same covariance structure; hence they provide more conserva¬ 
tive estimates of very extreme joint risks while maintaining asymptotic independence. Statistical 
inference is illustrated on extreme wind gusts in the Netherlands where a comparison to the Gaus¬ 
sian dependence model shows a better goodness-of-fit. In this application we fit the well-adapted 
Weibull distribution, closely related to the Laplace distribution, as univariate tail model. 

Keywords: spatial extremes; threshold exceedances; asymptotic independence; elliptical distribution; 
joint tail decay; wind speed. 


1 Introduction 


Extreme value analysis provides a toolbox for mode ling and estimating extreme events in univariate 
multivariate, spatial and spatiotemporal processes ( Colesl 2001 : Beirlant et ah . 2004 Davison et al 


2 OI 2 II . Principal objectives of spatial extreme value analysis are the spatial prediction of extremes and 


the extrapolation of return levels and periods beyond the historically observed range of data. A major 
distinction of dependence types can be made between asymptotic dependence when limu|i pr(Fv 2 (X 2 ) > 
u I Fxi (Xi) > tt) > 0 for two random variables Xi ~ Fxi and X 2 ~ Rvs and asymptotic independence 
when the limit is 0, provided the limit exists. Asymptotic independence can arise in environmental and 
climatic data for space lags or time lags when the most extreme events become more and more isolated 
in time, space or space-time. For many processes like wind gust speed or heavy rainfal l such behavior 
seem s plausible owing to physical limits, and it is corroborated by empirical findings (jPavison et al. 


I2OI3I ; Thibaud et aL . fic^ Ooitz et al. . 2015 1. In this paper, our objective is to construct asymptot¬ 


ically independent spatial processes that are flexible and tractable models with useful properties for 
modeling threshold exceedances. 

Models for asymptotically independent extremes must adequately capture the j oint tail decay rates 
in multivaria te distributions. A first in-depth ana lysis of joint tail decay was g iven bv iLedford and Tawn 
( 19961 1997 1. Closely related bivariate models ( Ramos and Ledford . 2009ll provide flexibility in the 
joint tail, yet an explicit definition of the probability density cannot be given when only one compo¬ 
nent is extreme and the generalization to higher dimensions suffers from the c urse of dimensionality. 
A mor e flexible characterization of multivariate tail behavior was developed by IWadsworth and Tawn 
( 2 OI 3 II . Useful models pertaining to this framework are obtained by inverting max-stable processes 
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( Wadsworth and Tawnl . l2012h . allowing for composite likelihood inference. Another a pproach that 
spans both asymptotically dependent and asymptotically independent data is presented bv lWadsworth et al. 
( 2014h . who model bivariate tails by assuming independence among the radial and the angular variable 
in a pseudo-polar representation. 

For lack of a unified theory of asymptotic independence, a variety of modeling approaches have 
proven useful in practice. They usually suffer from at least one of the following restrictions: joint tail 
decay rates are difficult to characterize; standard inference methods like classical likelihood are not 
available; only bivariate models are tractable and useful; the generalization to the infinite-dimensional, 
spatial setup is not possible; the univariate tail models prescribed by extreme value theory do not 
directly appear as marginal distributions in the model, necessitating marginal pretransformations that 
are not natural in the extreme value context. 

In the following, we present the novel Laplace model for multivariate and spatial extremes. It pro¬ 
vides a good compromise with respect to the aforementioned potential shortcomings. It is parametrized 
by a covariance function and closely related to standard Gaussian processes based on embedding an 
exponentially distributed variable for the variance. The resulting univariate distributions are of the 
Laplace type, and the univariate tails correspond to valid generalized Pareto limits of threshold ex¬ 
ceedances. Classical likelihood inference using threshold exceedances is straightforward. Joint tail 
decay rates and conditional distributions can be characterized in various ways. We use the terminol¬ 
ogy of spatial processes for simplicity’s sake, but an extension to the spatiotemporal context is possible 
through spatiotemporal covariance functions. The notion of a multivariate threshold exceedance is 
not uniquely dehned; here we will concentrate on four sensible choices for extreme value analysis: 
exceedances observed at one fixed site s, exceedances of a linear combination of values at D fixed sites, 
or exceedances of either the maximum or minimum value over D fixed sites. 

Section 2 gives a short exposition of some aspects of classical extreme value theory that are nec¬ 
essary to understand univariate tail models, their link to standard Laplace marginal distributions and 
the notion of asymptotic independence. Section 3 treats definition and inference of the asymptoti¬ 
cally independent Laplace model for threshold exceedances, whose tail behavior is characterized and 
contrasted with the asymptotic dependence case. An application to modeling of spatial wind gust 
extremes in the Netherlands in Section 4 is followed by the concluding remarks in Section 5. 


2 Extreme value theory 


Presentations going b e yond this short account can be found in iResnickI (|l987t) : iBeirlant et al.l (|2004ll : 
de Haan and Ferreira ( 2006ll . 


2.1 Univariate limit distributions 

The fundamental limit relation of extreme value theory establishes a three-parameter limit distribution 
for adequately rescaled maxima of an independent and identically distributed sample i = 1 ,..., n: 
if normalizing sequences On > 0, bn exist such that 

max a~^{Xi — bn) ^ Z, n —>• oo, (1) 

with a nondegenerate random variable Z, then the distribution of Z is necessarily of the generalized 
extreme value type with cdf 

Fziz) =exp(^- [1 + ?^]+^^^) > 

with parameters for position for scale ct > 0 and for shape When ^ = 0, Fz{z) is defined as 
the corresponding limit exp (—exp[—z]) l(z > 0). The tail index ^ is crucial to distinguish between 
Pareto-type (polynomial) decay for ^ > 0, exponential decay for ^ = 0 to a finite or infinite essential 
supremum and inverse polynomial decay to a finite essential supremum for ^ < 0. By transforming Z to 
Z* = 0.5 (1 + £,[Z — we normalize to a Frechet distribution Fz*{z) = exp(—l/[2z]) l(o,oo)(2). 
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Equivalent formulations in terms of threshold exceedances have proven useful both from a concep¬ 
tual and practical perspective. By absorbing the values of a„ and into the parameters a and /i, it 
is then useful to anchor practical modeling in the tail approximation 

pr(X > a;) (l = _logEz(x), (2) 

valid for large x. When C = 0, this tail probability is pr(X > a;) Ri exp(—[a; — m]/®"): corresponding to 
a generalized exponential tail. By conditioning on an exceedance over a high threshold u with u > fj, 
if ^ > 0 and u < otherwise, we get 

pr(X <a;|X>M)R;l — fl-l- J j = ct + ^{u — x> u, (3) 

where the right-hand side defines the cdf of the limiting generalized Pareto distribution for threshold 
exceedances x — u. 


2.2 Normalized marginal distributions 


When we treat the dependence structure in a multivariate setting, it is convenient to normalize all 
univariate marginal distributions to a parameter-free target distribution that appears naturally in 
the extreme value context. We will focus on standardizations that lead to limit distributions in 
H]) of the Frechet type Fz*{z) = exp(—l/[2z]) {z > 0) with tail index 1 and scale 0.5; we call this 
standardization the ^-transformation in the following. Here we slightly modify the usual normalization 
found in the literature (with scale 1) in order to establish the link to the Laplace distribution more 
easily in the following. In the literature, we often find X* = 1/[2(1 — Fjf(X))], leading to a Pareto 
target distribution with shape parameter 1 and lower bound 0.5 when Fx is continuous. Similarly, 
X* = —l/[21og(Fjf (X))] establishes the Frechet target distribution Fx*{x) = Fz*{x) when Fx is 
continuous. It is possible to use more general probability integral transforms for * , as long as they 
lead to nonnegativity and standard Pareto tail behavior, i.e., pr(X* > 0) = I and 2x x pr(W* > a;) —>• 1 
for X ^ oo. 

Concerning the behavior in the bulk of the target distribution, there is no other constraint than 
nonnegativity since the normalization of X* by a„ = n drives all non-extreme values from the bulk to 
0 in the limit. However, when modeling and estimating multivariate or spatial extremal behavior from 
a finite sample o f data. , we often also need to specify t he bulk behavior in the ta. r get di stribution, see 
Coles and TawnI ( 199 il l : IWadsworth and Tawn ( 2012 1: de Carvalho and Davison ( 2014h . for instance. 
For a useful choice of target distribution, we here focus on the following two properties: (1) a decreasing 
density on (0, oo): fx*{xi) > fx*{x 2 ) for any 0 < a;i < 3:2 < oo (2) exact Pareto tail behavior: 
pr(X* > x) = l/(2a;) for a; > u with u > 0.5 chosen as small as possible. These conditions yield a 
^-transformation ensuring a distribution of X* that is as close as possible to the asymptotic univariate 
limit since the convergence of maxima to Z* is equivalent to the convergence of the renormalized sample 
{X* /n, i = 1,..., n} to a Poisson limit process on [0, 00 ) with decreasing intensity da;/(2a;^) associated 
to the intensity measure A([a;, 00 )) = l/(2a:). The only distribution to satisfy these requirements is 
characterized by the density 

r l/(2a;^) if a; > 1, 

fx,{x) = lo.5 if0<x<l, (4) 

[0 otherwise. 

This target distribution is a mixture with probability 0.5 of a uniform distribution on (0,1) and a 
standard Pareto distribution on [1, 00 ). Its importance for the following developments com es from con¬ 


sider ing the distribution of log X*, whose density defines the standard Laplace distribution (|Kotz et al 

[2^, 

fiogX*{x) = 0.5exp(-|a;|). 

Therefore, the Laplace distribution can be considered as a natural standard marginal distribution for 
the modeling of multivariate extreme values. 
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2.3 Multivariate and spatial limits 

Univariate limits are readily generalized to the multivariate and spatial domain. For a sequence of 
independent and identically distributed copies {Xi(s)} {i = l,2,...,n) of a spatial process X = 
{X(s)}, we need normalizing sequences an{s) > 0 , bn{s) such that the convergence of componentwise 
maxima holds in terms of finite-dimensional distributions: 

max an{s)~^{Xi{s) — bn{s)) ^ Z{s), n —>■ oo, (5) 

with nondegenerate marginal distributions in the limiting max-stable process Z = {Z{s)}; we say that 
X is in the maximum domain of attraction of Z . The univariate convergence of marginal distribu¬ 
tions corresponds to the theory in Section EH We say that a vector X(s) = (X(si),... ,X{sd)) is 
asymptotically independent if the components Z{sj), j = are mutually independent. The 

limiting dependence structure of a multivariate vector can be characterized for ^-transformed margins. 
Convergence ([5|) is equivalent to the convergence of all univariate distributions in the sense of m 
together with the following condition: a (—l)-homogeneous measure rj with T]{tA) = t ^r]{A) exists 
such that, for any set A relatively compact in [ 0 , oo] \ { 0 } C K with ri{dA) = 0, we have 

tpT{X*{s)/t G A) ^ r]{A), t ^ oo, ( 6 ) 

for the marginally normalized vector X*(s). Loosely speaking, the vague convergence property ([5|) 
means that we can use the approximation pr(X*(s) G A) ^ foi' “extreme” sets A bounded far 
from the origin 0 . In practice, the homogeneity of rj provides a simple formula for the extrapolation 
of extreme event probabilities to more extreme yet hitherto unobserved events: pr(X*(s) € tA) « 
t“^pr(X*(s) G A), for t > 1 and an extreme event A. By transforming W*(s) to logX*(s), this 
relation can equivalently be written as pr(log X*(s) — t G A) exp(—t) fi{A) with 77 = 770 exp, t > 0 
and A C an extreme event. Finally, we point out that ? 7 (( 0 , 00 )) = 0 when X{s) is asymptotically 
independent, i.e., the limit measure 77 has a singular structure with all of its mass concentrated on the 
Euclidean axes. 


2.4 Tail dependence coefficients 

Bivariate tail dependence coefficients are useful summaries of tail dependence, and in a spatial context 
we can study their evolution with respect to the spatial lag between two sites. The tail correlation 
coefficient A of a bivariate random vector {Xi,X 2 ) is defined as 

A = limpr(Fx 2 (W 2 ) > u \ Fx^iXi) > u) 

u-\l 

= limpr(Fx2(W2) > u,Fx^{Xi) > m)/(1 - u) G [0,1]. (7) 


It is well-defined and symmetric in Xi and X 2 if {Xi,X 2 ) is in a bivariate maximum domain of 
attraction. Its value is 0 if and only if Xi and X 2 are asymptotically in dependent. If this is the 
case, more informat ion is provided by the residual coefficient p G [0,1] ([Ledford and Tawnl [1996; 


Draisma et al.l . 12004!) , where 


pr(-^r > X, X 2 > x) = £(x)x 


( 8 ) 


for large x with £(x) > 0 a slowly varying function such that i{tx)l£{t) —1 for a: > 0 and t ^ 00 . 
When p is positive, it corresponds to the tail index ^ of min(W[',X 2 ). For independent variables Xi 
and X 2 , we have p = 0.5. We remark that p is always equal to 1 when Xi and X 2 are asymptotically 
dependent, such that it is of interest to study either A or p, depending on the presence of asymptotic 
dependence. For the bivariate Gaussian distribution, the residual coefficient is (1 -f piin)/2, with pun 
the linear correlation coefficient, and covers the full range [0,1] of theoretically possible values when 
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Plin varies from —1 to 1. In the following, we exclude the case p = 0 corresponding to strong negative 
assocation of Xi and X 2 . The residual coefficient p can also be defined more generally for D > 2 by 

pr(^r > X,..., X^ > x) = (9) 

for large x if an adequately defined slowly varying function i exists. 


2.5 Prom modeling to inference of extremal dependence 


Since there is no natural ordering of multivariate data, different notions of threshold exceedances can 
be of interest with respect to the application context but also with respect to tractability of statistical 
inference. We focus on four useful choices that are common in the literature. Given a fixed threshold 


u, a vector a; is a marginal exceedance in Xi if xi > u, it is a sum exceedance if 

is a max exceedance if maxj=i^,,,^£) > u, and it is a min exceedance if miuj^i jjXj > u. For 


the corresponding exceedance sets, we write Ax^{u) = {x 


i(w) = {x I > u}, 


Xi ^ n}, ^sumV*^y — 1*^ I 

Araax{u) = {x \ maxj Xj > u} and Amin (it) = {x \ miTijXj > u}. The limit (jH) holds for these sets 
when u > 0. 

Whereas models for asymptotically dependent data have a strong theoretical foundation, there is no 
natural model class for describing how asymptotically independent data approach their limit. Using the 
limiting asymptotic independence would usually be too simplistic (by imposing independence among 
maxima of components, for instance), and it would considerably underestimate joint risks at the 
subasymptotic level. It is therefore sensible to focus on extremal dependence models that adequately 
capture the rate of convergence towards asymptotic independence. A simple approach consists in 
combining an appropriate model for univariate t ails with a Gaussian dependence mode l. Known as 
Gauss ian copula modeling ( Davison et al.l. 20I2f). this approac h has been explored by Bortot et al 


1 200fil for oceanographic data and by Renard and Lang (2007) for hydrological data, among others. 
Its extension to the spatial domain is straightforward. 

For exploratory analysis, it is useful to investigate estimates of tail dependence coefficients based 
on the empirical counterparts of © and (|S]). In the spatial context, one can study their variation with 
distance. For Gaussian copula modeling, likelihood estimators piin of the linear correlation pun, based 
on threshold exceedances, are readily defined and yield an estimate (I+piin)/2 of the residual coefficient. 
Likelihood estimation based on threshold exceedances is a standard approach for both asymptotically 
dependent and asymptotically independent data. After fixing an exceedance region A, we assume a 
parametric distribution Fa according to a parameter vector 6 for exceedances X \ X £ A Fa- 
The choice of threshold it in A is rarely straightforward since it is usually a trade-off between bias 
and variance. If a lower threshold decreases estimation variance by using more information from data, 
one also has to take into account the rate of convergence towards the asymptotic regimes used for 
modeling. The likelihood L for jointly estimating the exceedance probability pr(JA £ A) and the 
extremal behavior is based on the censoring of data x not falling into A: 


L{e;x) = 1ac{x) X (1 - pr(X G A)) -h lA(ic)pr(X G A)fA{x) 
where Ja is the density of the exceedance distribution Fa- 


( 10 ) 


3 Properties and modeling of exceedances of Laplace random 
fields 

3.1 Laplace random fields 

We denote W = {IU(s),s G K} a Gaussian random field defined on a nonempty domain K- Multi¬ 
variate distributions represent a special case by setting K = {!,... ,D} for D G N- We provide the 
constructive definition of Laplace random fields as a mixture of Gaussian fields based on embedding 
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an exponential variable with scale 2 for the v ariance of the Gau ssian field. This definition and basic 
properties can be fo und in the monogr aph of iKotz et al.l (1200111 . Multivariate Laplace distributions 
are also discussed in lEltoft et al.l ( 2006ll . 


Definition 1 (Laplace random field). For W a centered Gaussian random field and a random variable 
T > 0 with following an exponential distribution with scale 2 and independent of W, the random 
field 

X = {X(s), s&K} = {YW{s), s G K} (11) 

is called Laplace random field (or Laplace field in short). We call X a standard Laplace field if 
W is a standard Gaussian field. If W has a eovariance function or covariance matrix E, we write 
X ~ £(£). 

The distribution of Y is known as Rayleigh distribution with parameter 1. For standard Laplace 
vectors, we use the notation S* to indicate that E* is a correlation matrix. In the remainder of this 
paper, the default assumption is that E has full rank such that the inverse E“^ is well defined. The 
univariate marginal distributions Fx(s) of & standard Laplace field X are symmetric with mean 0 
and variance 2, and they have pdf fx{s){x) = 0.5exp(— |a;|) corresponding to the univariate Laplace 
distribution. The variable Y has pdf friu) = 2 /exp(—0.5?/^). The covariance function of X is 2E. If 
E(s) is the covariance matrix of lT(s) = (IF(si),..., 1F(sd)), the conditional distribution of X{s) 
given Y = y is Gaussian with covariance y^T,{s). Using the notation i/ = 0.5(2 — D), the pdf of X(s) 

is _ 

fx{s){x) = (2^)0.5n|s(^)|(i.5 (0.25a;'E(s)"^a;)°''''' (^y/x'Y{s)-^x^ . (12) 

Here, is the modified Bessel function of the third kind, see Equation (A.0.4) in iKotz et al.l (l200ll) . 
Whereas marginal densities are continuous but not differentiable at 0, the joint density fx{s) has a 
singularity at 0 for D > 1 such that fx(s){x) —>■ oo when ||a;|| —>• 0. Since our interest will be directed 
towards extreme events with ||a;|| strictly separated from 0, we can neglect this particularity. 

Laplace vectors are part of the larger class of elliptically contoured random vectors, which can be 
represented using a random radial variable R > 0 and a D x D deterministic matrix A that operates on 
a random vector [/, independent of R, with uniform distribution over the Euclidean unit sphere in M.^. 
Then RAU defines an ellipti c random vector w ith dispersion matrix E = AA'. For the multivariate 


Laplace distribution, we get (|Kotz et al.L 1200111 


fair) = 2i-0-25^r(0.5D)-V^-iiLo.5D-i(r), 


r > 0. 


For Laplace vectors, we get X{s) = YW{s) = YRy^r^g-^AU with Rw{s) and E(s) = AA' associated to 
the elliptical Gaussian vector lU(s). The dispersion matrices of the Gaussian vector and the Laplace 
vector are the same and the random vectors differ only by their radial variable. 


3.2 Asymptotic behavior of exceedances 

The univariate standard Laplace tails with pr(A(s) > x) = 0.5exp(—x) for x > 0 correspond to the 
univariate tail model (I2|) with ^ = 0, cr = 1 and p = — log 2. For exceedances above a fixed threshold 
u > 0, we get the conditional generalized Pareto distribution Tx(s)|x(s)>«(2;) = 1 — exp(—(x — u)) for 
X > u according to ©. 

We now characterize multivariate tail behavior of a Laplace vector denoted as X = (Ai ,... ^Xo) 
£(E) with E = Due to the elliptical structure of its distribution, weighted sums of 

the components of X are Laplace distributed: ~ C{u}'Yijj), and we have sum-stability 

with respect to the generalized Pareto tail family with ^ = 0, /r = — log 2 and = a/o/ScJ. The 
mean of the components of X follows a j 2 <D'^hh) distribution. When E = E*, then 

X* = exp(X) has standard log-Laplace margins ([4]) whose shape parameter is ^ = 1, and the geometric 
mean of components has tail behavior leading to a Pareto exceedance distribution with shape parameter 
( = D~^ Si<ji j 2 <D ^ 1 fo'' values above it = 1. This is in contrast to the multivariate standard 
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Pareto limits for asymptotically dependent distributions, where sums of components are again Pareto 
distributed with shape ^ = 1, see lFalk and Guilloul (|2008ll . 

A useful coefficient of tail dependence for Laplace vectors can be defined through the scale parameter 
cr = (X]i< 7 , 72 <-D ^ components, where are the entries of the correlation 


matrix S* associated to S in order to remove the effect of the marginal scale. 


When D = 2, this 


coefficient is a/ 2(1 + and varies between 0 (for a *2 = —1) and 2 (for (t *2 = !)■ When I? > 2, it 
is lower-bounded by 0 and reaches its upper bound D in the case of complete dependence = 1 

for all 1 < ji, j 2 < D). This coefficient allows us to see how pairwise coefficients relate to higher-order 
coefficients, which is not st raightforward for the stand ard coefficients A and p of extreme value theory 
defined in Secti on 12.41 see S chlathe r and TawnI (2003) for inequalities related to A in a multivariate 
context, and see lStrokorb et a'Lr ( 201lTll for spatial properties. 


Proposition 1 (Joint tail decay rates). The residual coefficient of (Ai,^ 2 ) ~ £(E*) with cr *2 = pun 
is p = 1/(1 -I- piin)/2. The joint tail decay is characterized by a conditional limit for x,y £ (0,1] when 
u 00 , 

pr(Ai > qi-^/u,X 2 > qi-y/u I Xi > A 2 > -)> [xyY/^'^P\ (13) 

where qp is the quantile of the standard Laplace distribution associated to the probability p £ (0,1). 
The multivariate residual coefficient defined in ^ is p = \J e'(S*)“^e for a Laplace random vector 
X ~ >C(E) with E* the correlation matrix associated to E and e = (1,..., 1)L 


Proof. For the bivariate results, we check the conditions stated in Theorem 2.1 of Hashorva ( 2010l) . 
which establishes joint tail decay rates for a large class of bivariate elliptical random vectors. The 
theorem states the following: if a function w{-) and the so-called Weibull tail coefficient 9 > 0 exist 
such that 

lim — exp(—t) Vt £ M, lim w{ct)/w{t) = Vc > 0, 

then p = ([1-1- piin]/2)®/^. In our case, we exploit that the Bessel function of the third kind has asymp¬ 
totic behavior K^{x) ^ l{2x) exp(—x) when x tends to infinity. Therefore, fnlr) ^ CRy/r ex.p{—r) 

with a constant cr > 0. Using I’Hopital’s rule, we obtain 

hm = hm = exp(-t). 


We have w{t) = 1 = with 9 = 1 yiel ding p = y^(l -|- piin)/2. The convergence (IT^ is a special 
case of Theorem 2.1 from Hashorv^ ( 20ir)[ l. where we exclude the cases a: > 1 or y > 1 to obtain our 
slightly simpler formulation in terms of conditiona l probabilities. The general multivariate residual 
coefficient is obtained in Example 1 of Nolde ( 2ni4h with a = 1 for the Laplace distribution. □ 


For given linear correlation coefficient pnn, the bivariate residual coefficient of Laplace distribu¬ 
tions is the square root of the Gaussian equivalent. The slower joint tail decay rate is due to the 
stochastic variance embedded in the Gaussian process. A bivariate Laplace vector with pun = 0 has 
uncorrelated components Xi and A 2 , but the corresponding residual coefficient l/-\/2 is different from 
1/2 corresponding to independent Xi and A 2 . The following proposition resumes the extrapolation of 
probabilities for exceedance sets. 

Proposition 2 (Extrapolation of exceedance probabilities). For X ^ ^(E) a D-dimensional Laplace 
vector, u > 0 a threshold and t = (t,... ,t) > 0 defining a translation X — t of the vector, we observe: 

pr(A -t£ Asnm{u)) = exp x t^ pr(A: G Asum(u)), 

pr(A - t £ A„:^(u)) = exp(-t/^/miJ)pr(X £ A„:j^(u)), 
pr(X - t £ Ainax(u)) ~ exp(-t/a)pr(X £ Amax(u)), 00, 

where a = Aymax(CTii, ..., (Jdd) in the last equation. When E = E* is a correlation matrix, we further 
have 

pr(X -t£ A^iniu)) ~ exp (^-y/e'{E*)-^e x t^ pr(A G Ai„in(u)), u -)■ 00 . 
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Proof. The exact tail decay rate for sum exceedances follows directly from the sum stability of the 
Laplace distribution, here with u) = e: 


( D 

pr(X - t € Asum(M)) = pr I ^ > Dt + u 

Vi=i 

= 0.5exp (^—{Dt + u)/\/e'Y,e^ 

= exp (^—Dt/Ve'Eej x pr(X S Asum(M))- 

Similarly, the exact rate of marginal exceedances is obtained by setting a; = (1, 0,..., 0). To derive 
the formula for maxima, we first consider the case where ^/ojf = ct > 0 is constant for j = 

Then X* = exp(X/cr) is of standard log-Laplace type with density (|3]) and satisfies the convergence 
By defining t* = exp{t/cr) and u* = exp{u/a), we find that both of the expressions u* pr(X* G 
^max('«*)) and rit*pr(X* G Amax(t*w*)) tend to 77(Amax(l)) = D as u* ^ oo, where 
is the extremal coefficient whose value is D under asymptotic independence. Since D > 0, we get 
pr(X* G ^max(t*M*)) ~ (t*)“^pr(X* G Ainax('w*))- In terms of X, t and u, we get the stated 
convergence rate pr(X — t G Ainax('u)) ~ exp {—t/a) pr(X G Amax('«))- In the case where the diagonal 
values of S are not constant such that = maxj-i^,,,^r> Ojj for jo G J with 1 < | J| < D, we start 
by observing that 


pr 


( maxXi„ > u 

joeJ 


< pr (X G ^max(w)) < pr 


( maxXi„ > u 
Jo&J 


^pr(Xj > u). 

jiJ 


(14) 


Since 0 < pr(Xj > M)/pr(maxj(,gj Xj^ > u) < pr(Xj > it)/pr(X^j > u) —?► 0 as u —?► oo for any 
ji G J and j ^ J, dividing the inequality (TT^ by pr(maxj(,g j Xj^ > u) shows that pr(X G vlniax(n)) ~ 
pr(maxjggj Xja > u) when u —>■ oo, and the convergence rate of maxima of X is the same as the 
convergence rate of maxima of the subvector Xj. Applying the result for constant diagonal elements 
(Jjgjg, jo G J to Xj then proves the statement for maxima. The asymptotic rate for min exceedances 
is directly linked to the definition of the multivariate residual coefficient p whose value is given in 
Proposition [H □ 


When S = E* such that margins of X are standard Laplace, then ([6|) holds for X* = exp(X) 
and the limit measure ij is zero for the sets exp(Asum('w)) and exp(Amin(u)) that lie within (0, oo). 
The decay rates are faster for these exceedance types than for and A^ax where decay rates are 
determined by the univariate marginal decay rate. Propositionshows that tpr(X*/f G A) tends to 
a finite positive limit when choosing i = i'/e's*e/D f Amin- 

The sum stability of elliptical distributions further provides exact formulas for related exceedance 
probabilities: pr(X G Asum(u)) = 0.5exp(—u/Ve'Ee) and A^^iu) = 0.5exp(—u/ ^au). To avoid 
plain Monte-Carlo approaches for calculating probabilities pr(X ^ x) = pr(Xj > Xj for at least one j) 
which can be too inaccurate or too slow in certain cases, we can instead exploit tailor-made algorithms 
for calculating multivariate normal probabilities <I>i;(a;) ( Genz and BretzL l2009ll . Univariate numerical 
integration then yields more accurate values based on 


pr(X ^x) = 1- / $E(a;/y)/v(y) dy. 


(15) 


3.3 Conditional distribntions 

The following proposition provides a basis for conditional simulation of Laplace random vectors con¬ 
ditional to the value of a subvector. When X = (Xi,X 2 ) U(E) is a D-dimensional Laplace vector 
with Xi = (Xi,..., Xd) for d G {1,..., D — 1}, we write E^-, 1 < j < 2 for the corresponding blocks 
of the dispersion matrix E. 












Proposition 3 (Conditional distributions). 

1. When X = YW ~ >C(S) with mixing variable Y > 0, Y^ ~ Exp with scale 2, conditioning on 
X = X yields the conditional density of Y given as 

fY\x=x{y) = c{x, E) X exp (-0.5[y^ + x'Y-^x/y'^]) , 

with the constant c{x, E) = (27r)“®'^^|E|“®'^//x(a;) > 0 depending only on x and E. For D = 1 
and E = (1), we obtain the conditional density 

fY\x=x{y) = 2° ®7r"° ® exp (-0.5[y2 + x^/y'^] + x) (16) 

whose mode is y/x. 

2. The conditional distribution of X 2 \ Xi = Xi is elliptical with density 

/x2|Xi=^i(®) = Co ((a; - fi)'t-^{x - /i) + ci) , (17) 


where, using v = 0.5(2 — D), 


5(c) 

Co 
T 

and Sn denotes the surface area of the unit ball in R" with si = 1 and Sn = 27r*’'®"'/r(0.5n), 
n > 1. The radius R in the pseudo-polar representation {X 2 \ Xi = Xi) = jl RAU with 
E = AA' has density 

fftir) = SD-d Cq + Ci), r > 0. (18) 


= (2l^(0-25v^)°'""A^(v^), 

pCO 

= SD-d 5(r^ + ci)r'°“'^“Mr, ci 

Jo 

= '^ 21 '^iiXi, E = E 22 — E2iEj^^^Ei2, 


— ^2^22 ®2, 


Proof. 1.) The change of variables {y,w) {y,x) = {y,yw) in fY,w = fY x fw yields 


fY,x{y,x) = y 




exp(—0.5?/^) X (27r) 


-0.5D|v|-0.5 


exp (—0.5[a:'E ^£c]/y^) . 


Using the formula /y|x=a:(y) = fY,x{y,x)/fx{x) leads to the given expressions. 

2.) The formulas for conditi onal elliptical dis t ributi ons follow from the related standard theory, see 
Corollary 5 and Section 4 of ICambanis et ahl (ll98T ). with g{-) given through the relation fx{x) = 
\Y\-°-^glx'Y-^x). □ 

One concludes from the density expression (fT6l) in Proposition [3] that the random variable Y \ 
A = X is concentrated around its mode i/x- Indeed, the ratio fY\x=x{x°'^~^^)/fY\x=x{x°'^) = 
exp(—0.5(x^+^'^ + x^~^'^) + x) tends to 0 exponentially fast for e > 0 when x —>■ 00 . Therefore, 
the distribution of the renormalized variable Yjx \ X = x converges to a point mass in 0 as x tends 
to infinity, and for X being large both Y and W must be large simultaneously. This is different from 
asymptotically dependent normal scale mixtures YW characterized by a mixture variab l e Y wi th tail 
index ^ > 0, where the most extreme events are due only to high values of Y, see lOoita (|2013ari . 

For simulating X 2 \ Xi = Xi, we could either directly calculate and simulate the elements of the 
pseudo-polar representation fi RAU, or we can exploit the latent Gaussian structure as follows. We 
first sample the mixture variable y according to the density /y|Xi=cci with a standard approach like 
inverse transform sampling based on the numerical integration of fY\Xi=xi- Given y, we sample W 2 


according to the conditional Gaussian distribution W 2 
covariance matrix E. 


Wi = Xi/y with mean E 2 iEj^^^a;i/y and 
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3.4 Likelihood inference 

We assume that the data process {X(s), s £ K} has been observed at D sites s = (si,. .., sd). For 
simplicity’s sake, we further assume that the marginal distributions of the observed random vector 
X(s) are of standard Laplace type. In practice, this requires marginal pretransformation according to 
an adequate univariate model ^x(s)- We adopt the exceedance-based approach characterized through 
the censored likelihood (nni). We have to fix an exceedance set A C to characterize the region of 
extreme events. Since we are assuming extremal dependence according to the Laplace random field 
model, we assume that the density of X(s) is equal to the multivariate Laplace density (IT^ for 

X G A. Denote pyi = pr(X(s) £ A) the exceedance probability for a Laplace vector. Given a correlation 
model with (parametric) correlation matrix S*(s), the likelihood contribution of an observed vector 
x{s) = {x(si),.. . ,x(s_d)) is 

L{E*{s);x{s)) = l(a;(s) ^ A) x (1 -pa(S*(s))) -f l(a;(s) £ A) x fx{s)ixis)). (19) 

When no tailor-made method exists for calculating the exceedance probability p^(S*(s)), a simple 
Monte-Carlo procedure can be applied by generating an independent and identically distributed sample 
of X{s) and by using the observed proportion of realizations inside the exceedance region A. We 
recommend a sample size of at least 10000, leading to a standard approximation error of around 
y^^/VTOOOO = y^^/100 for small pA- 


4 Application to wind gusts 


We illustrate modeling on daily maximum wind gust data from the Netherlands collected from 14/11/1999 
to 13/11/2008 for 30 meteorological stations, available for download from the Royal Netherlands Me¬ 
teorological In stitute (www.knmi.nl ) . Modeling extreme wind gusts is importa nt for applications like 
insurance risk dB rodin a nd Rootziil l2009i Mornet et al. ^ 2015 ), forest damage dPontailler et al.l. W97 ; 
IPhbtel . I 2 OO 5 I : iNagel et al.l . 2006 1 or wind farmin g (Seeuro and Lambert . 200^ Steink ohl et ah . 2^013 1. 


Recent studies on similar data (jEngelke et al 


[201a 


Einmahl et al 


2015t Oesting et al.L 12015 1 are 
based on max-stable models without challenging the assumption of asymptotic dependence. For the 
present data, we will show that the asymptotically independent Laplace model satisfactorily accomo¬ 
dates a situation that is inherent to asymptotically independent data: tail correlation estimates A tend 
towards lower values when the tail fraction used for estimation is decreased. Moreover, when using 
asymptotically dependent models, such data behavior makes the choice of thresholds for inferential 
methods difficult. O ther studies have shown that asymptoti cally independent models are preferable 
for wind speed data ( Ledford and Tawn . 1996t Ooitz . 2013bl) . hence we limit our analysis to asymp¬ 
totic independence models by comparing models based on marginally transformed Laplace fields or 
Gaussian fields. We removed a small number of observation vectors with missing components and 
retain observations for 3241 days. We conducted a preliminary analysis that showed moderate to weak 
day-to-day dependence of extreme wind gusts. Here we focus only on spatial modeling of intra-day 
dependence and neglect seasonal variations and clustering of extremes across several days. 

For estimating the parameters of the dependence structure with the likelihood (fTOl) . we use max 
exceedance sets Hniax('w). This yields an estimation strategy that is coherent for marginal modeling 
(using observations above the marginal threshold u) and dependence modeling (using observations 
within 4.i„ax('w) where at least one marginal threshold is exceeded). Since the parameter vector com¬ 
prising marginal parameters and dependence parameters can be of dimension up to 7 in our models, 
the joint estimation of marginal parameters and dependence parameters through numerical maximiza¬ 
tion of the likelihood is prone to numerical instabilities. We therefore estimate separately marginal 
parameters (through the independence likelihood) and dependence parameters. Since fitted univariate 
distributions often model imperfectly the actual data distribution, which is not unusual for spatial data 
due to their high dimensionality, we prefer to use the empirical probability integral transform for trans¬ 
forming data to univariate standard Laplace or standard Gaussian margins for estimating dependence 
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param eters. Computations have been carried out with the R programming language ( R. Core Teaml 


4.1 Data exploration 


Figure [T] shows the spatial setup of stations and information about sitewise quantiles for probabilities 
0.95 and 0.99. One can conjecture that wind gusts are slightly more extreme when sites are closer to the 
coastline; this conclusion can be drawn for both of the investigated probability values and is illustrated 
by the color code in Figure[T]highlighting the 10 highest and lowest quantile values among the 30 sites. 
Otherwise, quantiles seem stationary in space. Moreover, the ratio (go. 99 (s)/qo. 99 )/( 9 o. 95 (s)/qQ g^) 
with Qp the mean of the p-quantiles over all sites s is approximately constant, as can be seen from 
the information in Figure [1] where symbols + and x are of approximately the same size at each site. 
This justifies modeling the spatial nonstationarity of univariate tails through a nonstationary scale 
parameter in the Weibull tail model detailed in the following. Owing to the country’s flat surface, it is 
reasonable to assume that tail distributions vary only with respect to distance to the sea. Due to the 
rugged and irregular shape of the Netherlands’ coastline, we use one site’s distance to the segments 
connecting the coastal sites as a proxy for its distance to the sea. If a new site lies between this curve 
and the actual coastline, we set its distance value to 0. 

The left-hand plot in Figure [5] shows empirical estimates of the tail correlation function A(h) 
with respect to the distance between two sites, smoothed with local regression techniques (LOESS), 
for threshold values u given as 0.9, 0.95, 0.98, 0.99, 0.995. Estimates decrease when u increases, 
meaning that tail correlation weakens when we go farther in the tail, which is typical for asymptotically 
independent data. Note that the decision for or against asymptotic independence based on a finite 
data sample can rarely be made with absolute certainty in practice. Although tail correlation estimates 
move closer to 0 only very slowly when u increases, goodness-of-fit checks in Section 14.31 will confirm 
that the asymptotically independent Laplace model can reproduce such behavior. The assumption of 


asymptotic independence is therefore plausible, and the right-hand di splay of Figure [2] shows pairwise 


estimates of the residual coefficient p, calculated as the Hill estimator dHiliIl9'^ : Draisma et ah . 2004h 


based on the largest fc = 40 order statistics of min(X*(s), X*(s -|- As)) for all pairs of sites. From 
the pins in Figure [2] indicating the orientation of the spatial lag of site pairs, it is difficult to conclude 
on the presence of geometric anisotropy; we will therefore use model selection tools to decide on this 
issue. 


4.2 Marginal modeling 


The vast literature on wind s peed modeling has identified the Weibull distribution as the most adequate 
cand idate for univariate tails ( Stevens and Smuldersl . [l9^ISeguro and Lambertl . l2000HAkdag and Dinler . 
2 OO 9 I 1 . Therefore, we will not fit the classical tail model (l2|) and the related generalized Pareto dis¬ 


tribution ([3]) for exceedances to univariate tails, but instead the Weibull distribution whose density 
is X !->• l[o_oo)(a:)7(5“^x'^“^ exp[—{x/( 5 }^] with parameter constraints ^ > 0 and 7 > 0. The Weibull 
distribution extends the classical tail model with shape ^ = 0 (limit of Gumbel type) and mean p = 0 
by introducing a supplementary shape parameter 7; the classical model then applies to transformed 
data x'^. We account for spatial variation by allowing the covariate “distance to the sea in km” to 
modify the value of the Weibull scale parameter 6 = exp((5o -I- (5i x covariate). For a fixed threshold 
It > 0, the contribution of an observation x to the independence likelihood of univariate tail parameters 
is 


l(x < m) X (1 — exp [—{x/5}'’']) -I- l(x > It) X 7(5 '^x'^ ^ exp [—{x/5}'^] 


( 20 ) 


In preliminary studies, we further tried to include a location parameter in a way similar to the classical 
tail model ©, t.e. replacing x hy x — p with p < u \n ( 1201 ) , but we could not detect any substantial 
improvement in the goodness-of-fit and the numerical maximization of the likelihood became unstable. 
After fixing u to the empirical 0.975-quantile calculated from all observations, estimates and bootstrap- 
based confidence intervals of level 0.95 are given as follows: 7 = 1.72(1.55,1.86), (5o = 2.44(2.29, 2.52) 
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Figure 1: Measurement sites for wind gusts in the Netherlands. Coastal stations are connected with 
a dashed and piecewise linear curve. Each site is marked with a +-symbol, whose size is proportional 
to the empirical 0.95-quantile divided by the mean of all such quantiles, and a x-symbol, whose size 
is proportional (with the same constant as for 0.95-quantiles) to the empirical 0.99-quantiles, divided 
by the mean of all such quantiles. In each case, the sites associated to the 10 highest quantile values 
are marked in red, and the sites associated to the 10 lowest quantile values are marked in blue. 




Figure 2: Exploratory plots. Left: LOESS-smoothed empirical tail correlation functions A with respect 
to intersite distance, estimated for different thresholds u. Right: Empirical estimates of the residual 
coefficient p with respect to intersite distance and local regression curve (LOESS); pins indicate the 
orientation of the spatial lag. 
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m/s and = —0.0021(—0.0026,—0.0019) m/s, where confidence intervals have been calculated from 
a block bootstrap sample (block size 30, sample size 100). The effect of distance to the sea 5i is 
significantly different from 0 and indicates that the strength of wind gusts weakens when we move 
away from the coastline. We follow common practice in spatial extreme value modeling and use the 
empirical dis tribution below the th r eshold by assuming dense observat ions in the bulk region of the 
distribution ( Coles and Tawn . 1991 : Wadsworth and Tawnl . 2012 . 2013 1. 


4.3 Dependence modeling 

We consider Gaussian and Laplace dependence models with correlation functions of exponential, 
stable or Matern type. For the Matern model, we tried out a selection of regularity parameters 
V G {0.1, 0.15, 0.2, 0.25, 0.3, 0.4,1,1.5, 2.5, 5}. We further allow geometric anisotropy to accommodate 
a potentially different scale of dependence along one direction, for instance stronger dependence or¬ 
thogonal to the coastline for winds hitting land from the sea. By denoting 0 G [0,7r) a rotation angle 
and 6 > 1 a stretching along this direction, we replace the original bivariate column distance vector 
As by the matrix product 


fb 0\ /cos(0) 
1^0 l) l^sin(0) 


— sin(0)\ 
cos(0) j 


As = M As. 


Parameters 6 and b are estimated. We use Akaike’s information criterion AIC to select the best model. 
To make values of the censored likelihood (IT^ comparable for the Gaussian and the Laplace model, 
we first transform the original data x to x with uniform margins through the empirical transform for 
each of the 30 sites. If L{9;x) is the likelihood (1191) for either Laplace or Gaussian margins, we use 
instead the likelihood L in terms of x, 


L{e-x)=L{e-F-\x))x\{{f{F-\i,j)))-\ J = l,...,30, * = 1,..., 3241, (21) 

id 

where F is either the standard Laplace or standard Gaussian cdf, and / is the corresponding density. 

Table[T]sums up the estimation results, where the highest AIC values 2 log L(0; x)—2 dim(0) of both 
the Laplace and Gaussian Matern models were obtained for the shape parameter v = 0.25. As before, 
standard errors have been calculated through the block bootstrap approach. When calculating AIC, 
we consider Matern models with different values of v as different models, such that v is not counted 
for the model dimension. The shape parameters retained for the stable and the Matern class indicate 
nondifferentiable trajectories. Geometric anisotropy improves AIC, but the improvement is relatively 
small compared to the differences between the three classes of correlation functions. Since 6 > 1, 
we find that spatial dependence is weakest along the direction vector (cos(0),sin(0)) Ri (0.58,0.81), 
here calculated for the Matern Laplace model, whereas it is strongest when we move orthogonally 
to this direction. This direction approximately follows the orientation of the coastline and we can 
conclude that winds hitting land from the sea (or the other way round) are stronger dependent in 
space. Throughout, Laplace models outperform their Gaussian counterpart. Overall, we select the 
Laplace model with the anisotropic Matern correlation which obtains the highest AIC value. The 
Gaussian model with the same covariance family is second best. Since the isotropic Matern correlation 
function attains the value 0.1 approximately at distance x scale = x scale, the estimated 
dependence is strong across the study region. For an illustration of the goodness-of-fit, we look at the 
quantile-quantile plots of the distribution of the sum of the observation vector above the 0.95-quantiles 
in Figure m The sum distribution is a good choice since its tail decay behavior can be specified with 
an exact, nonasymptotic expression for elliptical distributions, see the results of Proposition [2] for the 
Laplace distribution. To facilitate the comparison of the displays for the Gaussian and the Laplace 
model, theoretical quantiles of the standard Laplace distribution are used and the data are transformed 
accordingly. For instance, the transformation of a data vector Xi = Xij, j = 1,... ,30 as defined in 
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Covariance 

Type 

e 

b 

scale 

shape 

log-L 

AIC 


G 


-(-) 

369.0(22.44) 

-(-) 

17420 

34830 


L 

-(-) 

-(-) 

376.6(18.43) 

-(-) 

17620 

35230 


G 

0.42(0.2) 

1.14(0.06) 

384.3(11.61) 

-(-) 

17430 

34850 


L 

1.02(0.15) 

1.09(0.05) 

381.9(7.7) 

-(-) 

17630 

35250 


G 

-(-) 

-(-) 

3294(5.14) 

0.53(0.01) 

17740 

35470 


L 

-(-) 

-(-) 

3299(47.49) 

0.50(0.01) 

17890 

35770 


G 

0.89(0.11) 

1.25(0.11) 

3302(3.97) 

0.54(0.02) 

17750 

35500 


L 

1.14(0.1) 

1.31(0.11) 

3299(58.8) 

0.52(0.01) 

17930 

35840 


G 

-(-) 

-(-) 

3149(381.4) 

0.25(-) 

17780 

35560 


L 

-(-) 

-(-) 

3087(297) 

0.25(-) 

17900 

35800 


G 

0.87(0.12) 

1.24(0.07) 

3384(136.4) 

0.25(-) 

18280 

36550 


L 

0.95(0.07) 

1.41(0.09) 

3698(121.7) 

0.25(-) 

18440 

36860 


Table 1: Estimated tail models for max exceedances of the Netherlands wind gust data. “Type” is either 
G for the Gaussian likelihood or L for the Laplace likelihood. Standard errors have been calculated 
through a block bootstrap approach. For the calculation of Akaike’s criterion (AIC; here values 
are rounded to 4 significant digits), the Matern shape parameter is not considered as an estimated 
parameter. 




Figure 3: Quantile-quantile plots for sum exceedances above the 0.95-quantile in the Gaussian (left 
display) and the Laplace (right display) tail dependence models with respect to the estimated Matern 
correlation models. In both cases, the theoretical quantiles correspond to a standard Laplace tail and 
the data quantiles are transformed accordingly. 


dl) into its corresponding standard Laplace quantile qi is done as follows for the Gaussian model: 


qi = F 


where F is the standard Laplace cdf and $ is the standard Gaussian cdf. In both cases, data points are 
aligned close to the diagonal, with no strong differences between the Gaussian and the Laplace model. 
Figure S] shows the difference of empirical and theoretical residual coefficients for these two models, 
with empiricial coefficients given as on the right-hand side of Figure [21 The local mean of differences 
has slightly positive values, which can at least partly be explained by second-order terms with respect 
to the asymptotic decay rate. Finally, Figure [5] contrasts the exact conditional exceedance probabilities 
Au(si, S 2 ) = pr(Fx(si)(Ar(si)) > u \ Fx(s 2 )(Al(s 2 )) > u) of the two models with the corresponding tail 
correlation estimates of data already presented on the left in Figure [2] Whereas the Laplace model 
convincingly reproduces the observed behavior in data, the Gaussian model is strongly biased towards 
too low values. 
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Figure 4: Difference of pairwise empirical estimates and theoretical values for the residual coefficient 
in the Gaussian (left) and the Laplace (right) tail dependence models with respect to the estimated 
Matern correlation models. Distances have been corrected for anisotropy according to the estimated 
anisotropy matrices such that dist = ||MAs|| 2 . 





Figure 5: Comparison of conditional exceedance probabilities A„ between data and htted Laplace and 
Gaussian models. From left to right for threshold values u given by {0.9,0.95,0.99}. The curves for 
data correspond to the empirical smoothed tail correlation estimates presented on the left of Figure [2l 
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Figure 6: Conditional simulations. Left: conditional to highest daily mean observation, observed on 18 
January 2007. Right: conditional to 10000 year daily return period of one coastal site. Points indicate 
the conditioning sites and their values. 

4.4 Conditional simulations and return levels/periods 

To illustrate conditional simulation for the Laplace model according to the algorithm described at 
the end of Section 13.21 Figure [5] shows examples for two scenarios. In the first one, we condition on 
values at the 30 sites observed on 18 January 2007 during the Kyrill storm, whose mean 32.4 m/s is the 
highest over the observation period. In the second scenario, we condition on the estimated daily return 
level Xcond = 51.0 m/s of one of the coastal sites with respect to a long return period of 10000 years, 
i.e., 1/(365.25 x 10^) = exp(—(xcond/'^)'^)- The gradient of marginal scale with respect to distance to 
the sea is well perceptible in both simulations. 

Regarding the calculation of return periods and levels taking into account spatial dependence, 
we first consider the return period of an exceedance of the level 51 m/s at at least one point in 
the Netherlands territory. Based on a fine spatial discretization si,..., S 345 of the Netherlands and 
the corresponding correlation matrix S* of the Laplace model, we calculate sitewise standard-scale 
values x{sj) = 1 — exp (—j = 1, •.., 345 with xo = 51. Using formula (TTK]) for numerical 
integration, the return period in years is 

^ 365.25x(l-/“<I.g.((^’-U5i).-T^U£345))/y)/v(y)dy) = 

where -F is the univariate standard Laplace distribution. Therefore, considering the whole country 
instead of a single site decreases the return period from 10000 to 376 years. We can further determine 
the return level associated to a return period of 10° years (c = 0,1, 2, 3,4) for the maximum observed 
over the whole territory. We numerically calculate the value Xq for which the function xq >—>■ |RP(xo) — 
10°| attains its minimum 0, yielding values (in m/s) 37.3, 44.0, 49.9, 55.1, 60.0 for c = 0,1, 2, 3,4 
respectively. 


5 Discussion and perspectives 

Due to its generalized Pareto tails and the sum stability of its multivariate elliptical distributions, 
the Laplace tail model has advantageous properties for extreme value analysis while remaining close 
to the Gaussian processes used in classical geostatistics. For conditional and unconditional simulation 
in high dimensions, the well-studied Gaussian techniques are the main ingredient. Standard maximum 
likelihood inference is possible. In contrast to the Gaussian copula model, it is easier to interpret 
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in the extreme value context since the elliptical structure arises for Laplace marginal distributions 
which, unlike the univariate Gaussian, are more natural for extremes. Moreover, the standardization 
to marginal Laplace distributions yields an exponential tail decay rate that is relatively close to the 
actually observed tail decay rates in environmental processes like wind speeds or precipitations. 

The goodness-of-fit checks for wind gust data revealed a better fit of the Laplace model in terms 
of Akaike’s criterion and bivariate joint tail decay behavior. We did not take into account temporal 
dependence among extremes. A simple possibility for capturing clustering behavior in time would be to 
model the time series of a summary statistics like the sum of the marginally norm alized observations at 
the observed time steps with models akin to the exponential ARMA processes of Lawrance and Lewis 
( 198 ^ . 


The Gaussian scale mixture construction of Laplace models opens perspectives to more complex 
models. The latent exponential variance variables could be dependent over time or vary over space, 
which would then necessitate a Bayesian framework to handle latent variables. Finally, other than 
exponential variables V could be substituted for the variance of the Gaussian field. Such modeling 
would lead to higher flexibility in the extremal dependence structure with respect to tail decay rates, 
yet densities and cdfs are usually not available in closed form, and typically the resulting marginal 
distributions are not natural in the context of extremes. 
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